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ABSTRACT 

This  thesis  discusses  the  round-off  errors  incurred  in  solving 
Linear  Programming  problems  by  the  Revised  Simplex  Method,  using  the  Product 
Form  of  the  Inverse  -  the  most  usual  method  of  solution.   The  matrices 
involved  are  assumed  to  be  sparse  as  is  usually  the  case. 

Also  discussed  are  the  reinversion  of  the  basis  matrix  and  iterative 
refinement  of  the  solution  obtained. 
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1 •   INTRODUCTION 

The  history  of  linear  programming  goes  back  to  the  discovery  of  the 
simplex  algorithm  in  1947  by  George  E  Dantzig,  and  its  implementation  on 
digital  computers  in  the  mid  l95o»s.  Little,  however,  has  been  written  about 
the  round-off  errors  produced  during  the  computation,  the  most  notable  exception 
being  a  thesis  by  Bartels  [l] . 

A  linear  programming  problem  may  be  defined  as : 

minimize  a  linear  objective  function 

T 
z  =  c  X 

subject  to  the  constraints 

Ax  =  b,  x  >  0, 

where  A  is  an  m  x  n  real  matrix  and  x,  b  and  c  are  real  vectors  of  length  n, 
m  and  n,  respectively. 

Two  approaches  to  error  analysis  are  possible.   The  forward  error 
a^al7s±s   meth°d  assumes  that  a  solution  x  +  Sx  has  been  obtained  and  the  size 
of  8x  is  then  estimated.  Alternatively,  the  backward  error  analysis  approach 
assumes  that  the  solution  obtained  is  the  exact  solution  of  the  slightly 
different  problem: 

minimize:   z  +  Sz  =  (c  +  5c)Tx 

subject  to:   (A  +  5A)  x  =  b  +  5b,  x  >  0 

and  estimates  of  the  size  of  5A,  5b  and  5c  are  made. 

Even  the  most  cursory  glance  at  input  data  for  a  linear  programming 
problem  will  convince  the  reader  that  the  information  is  far  from  exact. 


For  this  reason  backward  error  analysis  is  a  much  more  realistic  approach  and 
is  the  one  which  will  be  pursued  in  this  paper.   If  a  bound  can  be  given  on 
the  sizes  of  5A,  5b  and  5c,  and  these  values  are  within  tolerances  determined 
by  inaccuracies  in  the  data,  then  the  problem  may  be  considered  solved.   This 
places  on  the  user  the  responsibility  of  specifying  the  accuracy  of  his  data. 

A  feature  of  nearly  all  linear  programming  problems  is  that  the 
matrix  A  is  sparse,  with  very  few,  often  less  than  one  percent,  of  its  elements 
differing  from  zero.   This  fact  has  been  used  almost  from  the  outset  to  reduce 
the  amount  of  computation  and  storage  required  simply  by  not  storing  the  zeros 
and  omitting  calculations  involving  them.   It  is  also  possible  to  exploit  this 
sparseness  in  the  error  analysis. 

Control  of  errors  in  linear  programming  has  traditionally  been  done 
by  "tolerances".  Quantities  that  are  "very  small"  are  set  to  zero  and  numbers 
greater  than  some  small  negative  constant  are  deemed  to  be  positive.   Such 
tolerances  are  usually  set  by  "the  system"  but  may  be  changed  by  the  user  at 
his  peril. 

Efficient  solution  of  problems  seems  to  require  these  tolerances, 
and  techniques  for  their  automatic  control  have  been  developed  (e.g., 
Clasen  [2]).  However,  the  final  solution  should  be  free  of  any  such  arti- 
ficialities and  this  is  easily  achieved  by  iteratively  refining  an  approximate 
solution  to  purify  it  of  its  tainted  past. 

Chapters  3  and  k   cover  the  simplex  method  and  its  backward  error 
analysis  respectively.   Chapters  5  and  6  cover  a  reinversion  routine  that 
concentrates  on  retaining  sparseness  and  providing  error  control.   If  the 
solution  is  not  good  enough  after  reinversion  the  refinement  techniques  of 
Chapter  7  may  be  used.  Finally,  Chapter  8  shows  how  slight  imperfections  in 
the  solution  may  be  removed  by  perturbing  the  original  problem. 


2.      NOTATION  AND  PRELIMINARY  RESULTS 

It   is   assumed  that  all  operations   are  performed  in  normalized 
floating-point  arithmetic.      Define 


fl   (E) 


to  "be  the  result  of  evaluating  the  expression  E.  Following  Wilkinson  [9], 
the  following  laws  in  floating  point  arithmetic  are  assumed  to  hold  on  a  nor- 
mal machine,  providing  overflow  and  underflow  do  not  occur: 


fl  (a  +  b)  =  a(l  +  &1)  +  b(l  +  5  ) 
fl  (a  -  b)  =  a(l  +  5  )  -  b(l  +  5^) 
fl  (a  X  b)  =  (a  x  b)  (l  +  8  ) 
fl  (a  /  b)  =  (a  /  b)  (l  +  5.) 


for  any  two  machine  numbers  a  and  b,  the  quantities  |&.|,  which  are  relative 
errors,  being  bounded  by  some  small  constant. 

Two  degrees  of  precision  are  utilized: 
normal  precision,  in  which 


sil  <  £i 


and  extended  precision  where 


i'  -   2 


No  assumptions  are  made  about  the  relative  sizes  of  61  and  £g.  Normal  pre- 
cision is  used  throughout,  except  for  iterative  refinement.  Note,  however, 
that  many  routines  use  "double  precision"  for  all  calculations. 

It  is  assumed  that  operations  involving  a  zero  operand  are  exact; 


i.e 


• ) 


fl  (o  +  a)  =  fl  (a  -  o)  =  a 
fl  (o  -  a)  =  -  a 
fl  (o  X  a)  --   fl  (a  Xo)  =  0 
fl  (o  /  a)  =  o     (afo) 

In  practice,  however,  careful  programming  should  be  used  to  avoid  doing  these 
operations  at  all,  as  they  are  redundant. 

Upper  case  letters  are  used  for  matrices,  lower  case  letters  are 
used  for  vectors  and  for  the  elements  of  both  matrices  and  vectors.  Where  a 
sequence  of  matrices  or  vectors  is  used,  superscripts  are  employed  to  dis- 
tinguish them  and  their  elements  are  similarly  labeled.  Matrices  and  vectors 
of  small  terms  introduced  to  account  for  rounding  errors  are  prefixed  by  5 
i.e.,  5B  is  a  perturbation  applied  to  the  matrix  B.  All  vectors  are  column 
vectors  unless  the  context  implies  otherwise. 

Throughout  the  analysis,  the  l^   vector  norm 

v  I   =  max    v. 


is  used  and 


A  I   =   max  I  [  Ax 


n 
=  max  E    I  a. 


J=l 


ij 


for  any  n  x  n  matrix  A.    Since  sparse  matrices  are  being  used,  it  is  conven- 
ient to  define  the  function  V  by 


V  (t)  = 


1  if  t  f  o 


o  if  t  =  o 


for  any  scalar  t.  Thus  the  number  of  nonzero  elements  in  a  vector 


v  =  (v1,  v2,  ...,  v  )  is 


Lemma 


n 

■Zn  V  (V.) 
L=l     V  2. 


If   |0  |  <  €       k  =  1,  2,  ..._,  i 


then 


i 

n 

k=l 


u+y 


-  1 


<  i  e' 


where 


e     = 


iL-feJ-1] 


and   i  <  n      (l) 


Proof : 


Let  p   =   J    (1+9,  ) 


k=l 


>   J    (1  -  €)   =  (1  -  e)1 
k=l 


P    <  (1  +  e)1 


(2) 


-l 


p-1  <  (1  -  e)~     -1 


(3) 


Since 


[(1  _  e)i/2  .  (1  _  e)"i/2J 


1  -  (1  -  e)1  <   (l  -  €)_i   -  l 


Hence,  from  (l), 


<  0 


1  -  p  <  (1  -  e)1  -  1 


w 


Now  (2)  and  (3)  together  yield 


-l 


p-i  I  <  (1  -  e)    -l 


< 


H1- 


For  small  e,  e'  is  only  slightly  larger  than  6.   It  is  often  convenient  to  use 
e  instead  of  e  in  error  hounds. 

In  the  evaluation  of  the  inner  product  of  two  dense  vectors,  it  is 
necessary  to  estimate  the  error  in  calculating 

using  floating-point  arithmetic.  Let 


Then 


•  -0 


s.  .  =  f 1  (s.  +  x.y.)  i  =  1  2       n 


i+1 


=  S.  (1  +  0.)  +  Xiy.(l+Ai): 


where    |e  |  <  e,   |a.  |  <  € 

J*  1    ~* 


and     e  =  either  e  or  e 


p  =  s 
*     n+1 


n 
=  E  x  y  (1+A  f       ft    (i+e .  ) 
i=l         x    k=i+l     * 


s 


Using  the  lemma, 


=  Z   x.y.  [  1  +  (n  -  i  +  2)  7±] 
i=l 


where     \y  A   <   e  ' 


Z   x.y.  (1  +  b±)         ,    |Bil  <  (n  +  1)  €' 

i=l 


Now  consider  the  case  where  x  and  y  are  sparse  vectors.  Let 


n 

V  (> 

i' 


Sc  =  Z   V 


Then 


=  min  (q^,  q^) 


n 
and  q'   =  Z  V  (x.y.) 

i=l 


q'  <  q 


If  X.,  Y.  are  x.,  y.,  respectively,  where  x.yi  is  the  j-th  nonzero  term  in 


J   J 
the  sum,  then 


n  I  Z  x.y.)    n/  Z  x.y} 


q' 


=  Z  X.Y.  (1  +  A.) 

where  |  A .  | <(  q'  +l)€f 

J 


n 

Z   x.y.  (1  +  6.) 

.  ,   11      i' 
i=l 


where  |5.   <  (q+l)e 
1  l   — 


In  fact,  5.  =  A.;  however,  use  of  the  less  stringent  bound  is  justified, 
since  q'  is  hard  to  estimate  but  q  is  easily  obtained. 

It  is  sometimes  convenient  to  consider  one  of  the  vectors  (say  y) 
to  be  dense,  in  this  case: 


=  n 


so 


where 


1  =   "x 


^i?i  Xiyv =  A x±yi  u + 5i) 


|5.|  <  (^  +  l)  e< 
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3-   THE  SIMPLEX  METHOD 

A  brief  outline  of  a  version  of  the  Revised  Simplex  algorithm  is 
given;  a  much  more  detailed  analysis  may  be  found  in  Dantzig  t3l  or  Hadley  [6] 
The  problem  to  be  solved  may  be  put  into  the  form: 
minimize 


T 
z  =  c  x 


subject  to 


Ax  =  b 

x  >  o. 

Let  B  be  the  basis  matrix  [P  ,  ...,  P  1  where  Py  ,  . . .,  Py  is  a 

1        m         1        m 
set  of  m  linearly  independent  columns  of  A.  Then  a  basic  solution  is  given  by: 


*  -1 

x  =  B  x  b, 


r  1T 

where  x  =  L  x  ,  •  •  ■ ,  x  j  , 

1        m 


since  the  inverse  exists;  this  solution  is  also  feasible  if  x  >  0.  The  suc- 
cess of  the  Simplex  method  depends  on  the  fact  that,  if  an  optimal  solution 
exists,  there  is  a  basic  feasible  solution  which  yields  the  optimum  value  of 
the  objective  function.   It  is  therefore  necessary  to  consider  only  the  basic 
feasible  solutions  in  the  search  for  optimality. 

An  algorithm  for  solving  the  problem  by  the  Revised  Simplex  Method 

is  as  follows : 

(a)   Find  an  initial  basic  feasible  solution,  or 
indicate  that  none  exists. 
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(t>)      Change  the  basis  until  optimality  is 
reached,  which  is   achieved  by  a  series 
of  simplex  iterations.      The  steps  re- 
quired for  each  iteration,    starting  with 
some  basis  B  =  [P ,  P ,    . . . ,   p ]    are : 
1         2  m 


Step  1:       Set  yT  =  [^    ,    cy   ,    ...,   cy   ] 

12  m 


Step  2:        Compute  7T      =  7T  B_1. 


T 
Step  3:        Compute  d.  =  c.   -tt      P. 

J    J       J 

for  j  =  1  (1)  n,  J  f  {vi,  v2,  ...,  vj 

Choose  an  index  s  such  that 

d  =  min  {d.  :  d.  <  0  } 
s        3  J 

If  no  d.  <  0,  the  solution  is  optimal. 

J 


Step  k:       Compute  y  =  B~  P 


Step  5 :   Find  an  index  r  such  that 


x     mm 
r 


yr 


in  I'M 


If  no  y   is  >  0,  then  the  solution  is  unbounded 
Choosing  the  minimum  of  x./y.  assures  that  the 
new  basic  solution  x'  will  be  feasible. 


Step  6:   Change  the  basis  to 


B'  =  tpv  >  "•>   pv   '  Ps>  Pv   '  '"'   Pv  ] 

1        r-1       r+1        m 
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Step  7:   Compute  the  new  solution,  x'  ,  using 


x"  =  (B*  )_1  b. 


Since  three  sets  of  equations  axe  solved  using  the  matrix  B,  it  is 
advantageous  to  have  B~J"  available.  If  B  "  is  known,  then  the  inverse  of  the 
new  basis  (B' )"  is  easily  obtained,  since 


(B')'1  =  (B'^BB'1 

=  (B"1  B')"1  B'1. 


B'  is  identical  to  B  except  that  the  r-th  column,  Py  ,  is  replaced  by  Pg,  so 

r 
B~  B*  takes  on  the  simple  form: 


B^B'  . 


o 


•    o 


o 


o 


'm 


-1 


where  y  =  [y1,  y2>  '  "  '   ym    ~  B   Ps' 
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Such  a  matrix  is  called  an  elementary  column  matrix,  or  simply  an  elementary 
matrix.   It  is  easily  inverted  to  give: 


(b"M  "  = 


-y2/yr 


V 


yv 


~yJyi 


say,  and 


(B1)"1  =  T  B_1 


Glearly,  if  the  starting  basis  is  B  ,  and  IT,  T  ,  . ..,  T5  are  the  elementary 
matrices  produced  during  p  iterations,  then  the  inverse  of  the  current  basis 
B  ,  say,  is  given  by: 


"1  =  T9   n^"1   ...   T2  T1  B  _1 


It  is  usual  to  take  B  =  B    =  I,  so  that 

o    o      ' 


_1  =  T*  T5"1  .  .  .  T2  T1. 


If  a  unit  matrix  cannot  be  found  from  the  columns  of  A,  then  one 
must  be  manufactured  using  "artificial"  columns.  When  the  corresponding 
artificial  variables  are  all  removed  (or  reduced  to  zero)  by  subsequent 
iterations,  a  basic  feasible  solution  has  been  found. 


11+ 

When  the  matrix  B  ~  is  not  calculated  explicitly,  but  only  Or  , 
P 

rpP"1   ,  # .   t   t1  are  stored,  the  method  is  called  the  product  form  of  the 
inverse.  We  will  restrict  our  analysis  to  this  method,  since  it  is  used  more 
extensively  than  any  other  method  in  large  scale  LP  algorithms. 

In  order  to  reduce  round  off  errors,  storage  requirements  and 
running  time,  it  is  desirable  from  time  to  time  to  start  afresh  and  recal- 
culate the  inverse  of  the  current  basis.  This  technique,  called  re inversion, 
is  discussed  in  Chapter  5. 
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k.      ERROR  ANALYSIS  OF  THE  PRODUCT  FORM  OF  THE  INVERSE  METHOD 

In  this  chapter  the  errors  introduced  during  one  iteration  of  the 
Simplex  method  are  considered. 

Suppose  that,  after  p  iterations,  the  approximate  inverse  of  the 
basis  B  is  given  in  product  form  by: 


TP  T5"1   ...   T2  T1 


where  each  TJ  is  an  elementary  matrix.  Define  an  error  matrix  E  by: 


PP  TP_1  ...   T2  T1  (B  +  E)  =  I. 


(1) 


Suppose  also  that,  in  performing  one  iteration,  the  r-th  column 
of  B  is  replaced  by  a  non-basis  column  of  A,  Aq  =  w  say,  giving  a  new  basis 
matrix  B'.  A  new  elementary  matrix  TP+1  is  formed  which  satisfies: 


TP+1  TP  TP-1  m    t  ^   T2  Tl  (B,  +  E,}  =  i      (2) 


The  round  off  errors  in  (l)  and  (2)  are  measured  by  ||e||  and  ||e'||, 
respectively.   The  problem  is  to  find  a  bound  for  | |e' | |  in  terms  of  | |e| | . 
To  simplify  the  notation  assume,  without  loss  of  generality,  that 
the  j-th  elementary  matrix,  T°,  is  a  unit  matrix  except  for  its  j-th  column 
and  that  in  the  j-th  iteration,  it  is  the  j-th  column  of  B  that  is  replaced. 
Thus,  pivoting  occurs  down  the  diagonal.  This  can  be  achieved  by  suitable 
interchange  of  rows  and  columns.  In  practice  the  actual  pivots  may  be  re- 
corded by  saving  the  pivot  column  number  corresponding  to  each  component  of 
the  basic  solution.   Thus  set  the  j-th  elementary  matrix  to  a  unit  matrix 
except  for  the  j-th  column  which  is 
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JlT 


mj,  ti^,  ...,  ti^] 


Since  T  is  obtained  by  inverting  an  elementary  matrix  whose  j-th  column  is 
y  say,  define  the  vector  T\      to  be  given  exactly  by 


i  =  l(l)m,  i  |=  j 


i  -  3 


Whenever  the  T  is  used  in  calculation,  the  error  introduced  by  the  divisions 

will  be  included  in  the  error  analysis. 

Either  of  two  different  methods  may  be  used  in  calculating  T]  ,  in 

each  of  which  t  =  -y.  and  y.  =  -1  are  set  first.   One  method  is  then  to  divide 

J      J 

each  element  of  yJ  by  t.   The  other  is  to  multiply  each  element  by  t   .   The 
latter  method  will  be  faster  if  division  takes  longer  than  multiplication; 
it  does,  however,  produce  two  round  off  errors  for  each  element.   It  is 
assumed  that  the  latter  method  is  used. 

One  step  in  the  Simplex  algorithm  is  to  compute 


y  =  B  "  w 


for  whi»-h  the  approximate  solution  is  given  by 


y  =  f 1  (TP  T9'1    .  .  .  T2  T1  w)  (3) 


Setting 


1 
a     =  w 


17 


and  calculating  successively 


k+1         k  k 

a*  x  =  f l  (tV)    (k  =  1,  2,  ...,   P)     (h) 


yields 


+1 


y  =  a? 


Equation  (k)    is  equivalent  to 


fl  (a  +TliO^)    ±=lf    2,    ...,  m 

™k+l  i 

fl^kak)         i  =  k       ^) 


Rearrangement  of  (5)  to  include  the  effect  of  computations  on  previous 
gives : 


i-1 
fl  (Z        o£    T]*)  ,     ±<£     (6) 

CC 

fl  K  +k£  ak  *n±>       >  ±  >  ' 


In  particular,  setting  1  =  £,   yields 


°£  ■  fl  <wi  +Sl  °{  Tlj>-  (7) 


For  convenience  in  notation,  let  p,  =  a  .  Then 


pi  "  fl  <wi  \£  pk  ^i)  (8) 


a 


a:.  : 
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r.     =   «? 
1  1 


+1 


fl  (Li     \    1?  '2P  (9) 

n (wi + L  \  <]i>* 


Now  let  L  be  the  matrix 


L  = 


-Tl 
41 


2 
1 

3   "'3 


<        1 


1     2 

'm    'm 


Then  (8)  may  he  written  in  the  form: 


L  p  =  w 


from  which  3  is  uniquely  given  by 


3  =  L   w. 


Using  the  relations  given  in  Lemma  1  of  the  Appendix,  it  can  be  seen  that  3 
exactly  satisfies 

(L+6L)3=w+Sw, 


where 


19 


where 


(I* '" 


l»yl       <      (^  +  3)^ 


«i         =  iz1  v  cn±) 
k=i        1 


le  w±|         <     q.    |WjL|   e£. 
If       A  w  =  L  3   -  w 
then  Aw=§w-6L.   p 


11**11  <  ||e  w|  |  +  ||5l||   ||p|| 

<qL    lk||^+    ||6L||    ||?||  (10) 


qL  =  max       q_. 


i  x 


If  WL   =  max  |T]?| 

1  <  j  <  i  <  m  1  (X1) 


then  ||b  L||    <  e£      qL      (qL    +  3)   ^  (12) 


and  hence 


MMI    <     e{   (qLNw||    +   qL   (qL  +  3)   ^IIpH  )  (13) 
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Now  estimate  the  error  involved  in  calculating  y  using  (9) 
Let  U  be  the  matrix 


U  = 


<    < 


•n: 


T1J 


Tl 


1  2 

Tl  Tl 

'p+1      'p+1 


1  2 

Tl  T) 

m  m 


Tl? 


T| 


T! 


p+1 


in 


'1 


provided  that  p  <  m.   Then 


y  =  fl  (U  7) 


where 


7  =  [3-l,  P2,  .-.,  Pp,  wp+1,  ...,  wJT 


(1*0 


By  using  the  results  of  Lemma  2  in  the  Appendix,  it  can  be  shown  that  y 
satisfies 

y  =  U  (7  +  6  7) 


where 


7  =  [6  Px,  6  p2,  ...,  6  B  ,  Aw _  ,  ...,  AwJT  say, 
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67  < 


X 


u 


-1 


1-2.J  q^+2)  ^  IIU"1!! 

[€i  <u<V2>  "a  l|yH  +el  IMfl 


p. 


qTT  =  max  / 


,  i  >  p 
,  i  <  p 


(15) 


and 


=  the  maximum  number  of  nonzero  elements  in  any- 
row  of  U. 


"u       PS   (luijl)# 


1>J 


If  6  P  =  (8  31,  6  P2,  ...,  6  0  ) 


and  A  wf  =  (A  w   ,  ...,  A  w  ) 


then 


8  3||  <  | |6  7\\    and  ||  A  w« | |  <  ||6  7\ 


The  perturbation  6  P  of  3  gives  rise  to  a  perturbation  L  6  P  in  w  ,  ..., 
Therefore,  combining  (13)  and  (15) 


w 


y  =  f 1  (TP  TP_1  ...  T2  T1  w) 


=TP  T5"1  ...  T2  T1   (w  +  v) 
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where  ||v||   <  max  {eJ_(qLl|w||    +  qL(qL+3)  mJ  |3|  l)+|  N  |  |  |  63|  U  I  67|  |) 

(16) 


and  L   < 


Let 


stt 


To  simplify  (l6),  let 


M  =  max   (|n*|,  1)  (17) 

i,k 

Then  ^  <  M,  M^  <  M.  (18) 


v     k 
q  =  max  ^  V  (T^) 


Then  |  |v|  |    <  e|   (q|  |w|  |    +   q  (q+3)   M   | |p| |  ) 

+  £{  qM  ||    U_1|| 


1    •   2  £[   q  (q+3)   M    |  |U 


—  [q   (q+3)   M    ||y||    +    MpII    +   q||w||]  (19) 


-li 


The  only  quantity  in  the  above  inequality  which  is  unknown  is  j  |u  |  \; 
fortunately  approximations  to  the  elements  of  U  "  are  found  during  the  cal- 
culation of  the  T\   vectors. 

Consider  now  the  calculation  of  the  T]  vectors  forming  the  approximate 
inverse  of  3,  and  consider  3  to  be  a  unit  matrix  except  for  the  first  p 
columns.   If 


23 


and 


then 


ce- 


ll 


a 


21 


a- 


a 


12 

1 
22 


or 


a 


lp 

l 
2p 


B   = 


a 


a. 


pl         "p2 


a. 


'PP 


a 


ml 


a 


m2 


a 


mp 


k+1 
a±3 


for  i  =     1(1)  m       ,        j     =     i(i)p 


i  f  k 
i  =  k 
k     =  l(l)j-l 


(20) 


fl   (J 


d-i 


:=i 


<  ii> 


a 


ij 


k5i 


fl  (a1. 


i<  J 


i>  j 


(21) 


and  the  T\   vectors  are  given  "by 


*nv 


fl  {l/o?..) 

33 


if 


1   t   J 


1   =  3 


(22) 


Let  r  be  the  matrix 


2k 


r 


:■■ 


11 


a 


a 


12 

2 
22 


1 
a13  • 


Q 


Q 


23 

3 
33 


p+1,1 


aJ+l,2 


a: 


'ml 


a 


m2 


a 


a 


ip 

2 

2p 


ccx 


PP 

p+1,  p  1 


a 


mp 


Then 


k=i     kj    'i 


(UT)..   =     < 


Zn  ak .  T]k  +  a1 
£1     k,n     i         i: 


ij 


i  <  P,   J  <  P 

i  >  P,    J  <  P 

i  =  J  >  P 

otherwise  (23) 


But,    for  i  <    j  <  p 


L  <i  iS  ■  °jj  Td  +  £  %  ^ 


and,    from  (21), 


0( 


id 


3-1 


■Si  \j  il  +  .Si  »'ij  %  T|I 
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where 


I    hi    I   <  <%  +  2)   Si 


-k  =kkv(ni) 


£«£    <    -«L    ^IA   +  £    £««    if.  (A, 


Similarly,      for  i  >  p, 


A  <j  *i  +  4  -  -h  +  «i,  ^  +  £  &  <£  ^ 


1 


+  A    .      Gf  (25) 


but,    from  (22), 


aiJ     +    4  *!       -    2  sij  aij  *  <  J  (26) 


^  c^    n^  =  i  +  a?.  n^  bj!. 

where 


Hence 


ij    '   -      1 


ur    =      i  +  f 


where  F  is  a  matrix  of  small  terms  caused  by  round  off  errors 


If    q  =  max   L  V  (T)*) 

1    K--»-        ! 


M  =  max     (  |T1*|,  1)  , 
1      x 


as  before,  G  is  the  matrix 
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G  = 


4l    °t 


a 


°13 


22    a23 


1 
°1T 


a, 


2p 


c^  1  OP'} 
p-l  p-l  p-l  p 


cP 
PP 


and  H  is  a  p  X  p  matrix  given  by 


H .  .  -  aV . 


Then 


k=i     kj   i   ij       ij   ij 


j-l 


'«-\£    < 


Tl.   6.  . 
J   i   ij 


i  =  J  <  P 


i  <  j  <  p 


+  2  a?.  s!  .  +  a.  .  A.  . 

i  >  P,  J  <  P 


otherwise 


Hence 


fti'Jsld&^^i'  K\) 
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:en'  .£,   la"?. I  +  e,'  f  la1 

1  .1^1   i.i 


'1  j^L   '  ij 


Thus 


and 


Now 


<  ^  [q  (q+2)  M   |  |g|  |  +  2  |  |h|  |  +  ||b||] 


Ml  <    e[   U  (q+2)  m  ||g||  +  2  ||h||  +  ||b||]     (27) 


||r||  <  max  {  ||G||,  ||B||  +  1} 


u"1  =  r  (.1  +  f)"1 


u'^k   llrl 


(I  +  F)'1  || 


<     j |r|  I    provided  ||f||<  1 


(28) 
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Substituting  the  bound  for  | |u~"    into  equation  (19)  gives  a  bound  for 


It  has  been  shown  that 


2  _1 


y  =  fl  (T9  T9'1    .  .  .  T  T1  w) 


=  t5  T9'1    .  •  •  T  T   (w  +  v) 


It  is  now  possible  to  estimate   |e'||,  where  E'  satisfies 


rP+1  T9    •  •  •  T2  T1  (Bf  +  E')  =  I. 


The  new  elementary  matrix 


T9       has  as  its  p  +  1  -st  column  the 


vector  T|    where 


f  "yi  /  Vi 
I    1  /Vi 


i  f"  P+l 
i  =  p+l 


Thus 


p  zeros 


)+l 


T^1  y  =  [  0,  .  .  .  ,  0,  1,  0,  .  .  •  ,  0] 


=  T9*1   TP  •  •  •  T     T1  (w  +  v) 


since 


w  is  the  p  +  1  -st  column  of  Bf,  v  is  the  p+l  -  st  column  of  E' 
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For  all  other  columns  B*  and  B  are  identical,  and,  since 

^  .  .  .   T2  T1  (B  +  E)  =  tP+1  I 

=  ^+1   , 

(B1  +  E')  and  (B  +  E)  are  the  same  except  for  column  p+1;  i.e.,  E'  is  E  with 
the  p+1  -st  column  replaced  by  v.  Hence 

Ne'II  <||e||    +    ||v|| 

Finding  a  bound  for  ||v||  is  not  as  difficult  as  equations  (19) 
and  (27)  would  seem  to  indicate  since  the  only  quantities  required  that  are 
not  recorded  are  the  elements  of  G.   (3  is  just  the  last  column  of  G) .   These 
values  are  obtained  during  the  calculation  of  y  and  may  be  used  to  update 
||g||  from  iteration  to  iteration. 

This  bound  depends  heavily  on: 

1.  The  sparseness  of  the  T)  vectors; 

2.  The  size  of  their  elements. 

It  is,  therefore,  very  important  to  reduce  the  magnitude  of  the  T)  elements  by 
controlling  the  size  of  'pivots'.   Consideration  also  should  be  given  to 
maintaining  spareseness. 
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5-   REINVERSION  OF  THE  BASIS  MATRIX 

During  the  execution  of  the  simplex  algorithm,  the  choice  of  pivot  is 
restricted  at  each  iteration  by  the  need  to  obtain  a  significant  improvement 
in  the  value  of  the  objective  function  and  to  maintain  feasibility.  Little 
can  be  done  to  control  round-off  errors  and  the  sparseness  of  the  inverse.   In 
recalculating  the  inverse  a  flexible  pivot  selection  agenda  can  be  adopted, 
and  the  errors  and  density  be  thereby  greatly  reduced. 

Before  discussing  the  general  reinversion  algorithm,  the  decomposi- 
tion of  a  sparse  matrix  A  by  two  different  methods  will  be  considered.   In 
both  cases  pivoting  occurs  down  the  diagonal. 
(l)  Product  form  decomposition. 


1  2 

A  =  EXE   . 


E11"1  En 


where 


ri 


EJ  = 


el3 
e2j 


1  e 


e  .  . 

33 

e  •  -i  • 


nj 
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(2)      Decomposing  A  into  a  product   of  lower  and  upper  triangular  matrices 


P 


11 


L  = 


i21         ^22 


'31 


^32        ^33 


(jnl  ^2       ^n3 


1 


nn 


and 


U  = 


\2  \3 

1  U23 


•    '      "l» 


u 


2n 


u 


h-l,n 

1 


so  that  A  =  LU. 


Having  decomposed  the  matrix,    inversion  is   easy,    s 


nice 


is6)'1 


-e    ./e .  . 


'    V< 


jo     i 


-e    ./e  .  . 


.      1 


where 


L"1   =  Ln  Ln"X  .  .  .  L2  L1 
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LJ  = 


1   • 
1 


-i  i/i,. 


and  similarly  for  U   .  Note  that  the  pattern  of  zeros  and  nonzeros  in  E  , 
LJ,  UJ  correspond  exactly  to  that  of  the  matrices  from  which  they  were 
derived. 

The  density  of  the  product  form  decomposition  is  at  least  as  great 
as  in  the  LU  form,  as  can  be  shown  thus : 


^   = 


JJ 


>1,J   1 


nj 


'id 


1   e.  .  , 


TJ  (say). 


112? 
Then  A  =  S  T  S  T 


sn-l  Tn-1  sn  ^ 


For  i  <  j   ,   T1  SJ  =  SJ  T1 


So  the  product  of  the  S's  and  T's  may  be  rearranged  to  give: 


A  =  S1  S2 


on-l   n  1   2 
S    S  T  T 


T^1  Tn. 


.1   „2 


S  '  S  '  *  '  •  >   s  are  all  lower  triangular  matrices,  so  L  =  S1  S2 
is  lower  triangular. 


1   ? 
Similarly  U  =  T  T 


n 


33 


.  .  T  is  upper  triangular,  and  LU  is  the 
(unique)  decomposition  into  lower  and  upper  triangular  parts. 


L  .  S1  S2 


'11 


'21   1 


'nl 


623 


'n2 


.1 


'n-l,n-l 


'n  n-1  1 


'I 


nn 


'11 


e     e 
21   e22 


e31   e32    e33 


e  -,        e  _    e 
nl    n2     n3 


.  e 


nn 


Hence  L  contains  the  same  number  of  nonzero  elements  as  in  S1,  S2, 
(excluding  unit  diagonals  which  need  not  be  stored). 


1  2 

U  =  T  T 


T 


n 


.'.  U"1.^11)-1  (T11"1)-1  .  .  .  (T2)-l  (Tl)-1 


.,    S 


3U 


-e 


-e 


In 
2n 


-e 


-e 


ln-1 
2n-l 


'12 


-e 


n-ln 


"U 


■ 

1 

"el2 

"el3 

• 

• 

• 

-em 

1 

"e23 

• 

• 

• 

"e2n 

1 

• 

• 

• 

~> 

• 

'1 

* 

n-l,n 

1 

- 

- 

-i  O  yi 

Therefore  the  number  of  nonzero  elements  in  T  ,  T  ,  .  .  . ,  T  is 
equal  to  that  in  u"  (again  excluding  the  unit  diagonals). 

Since  the  density  of  U~  will  he  greater  than  that  of  U,  due  to 
"fill  in"  during  the  inversion  procedure,  it  may  be  concluded  that  the  decom- 
position into  elementary  matrices  will,  in  general,  give  rise  to  more  non- 
zero elements  than  will  LU  decomposition.   (TltiS  ignores  the  possibility  of 
cancellation,  which  might  produce  zero  elements  from  sums  of  nonzero  ones; 
however,  this  should  not  often  happen  in  practice). 

It  is  therefore  desirable  to  use  an  LU  decomposition  method  for 
reinversion,  not  only  to  reduce  the  rounding  errors,  but  also  to  reduce  the 
number  of  T|  elements  representing  the  inverse. 

It  has  been  shown  that  a  triangular  matrix  decomposes  readily  into 
elementary  matrices.   It  is  therefore  desirable  to  sort  the  rows  and  columns 
of  the  basis  to  form  a  matrix  which  is  as  close  to  a  triangular  matrix  as 
possible.   The  sorted  basis  is  then  partitioned  as  follows: 
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M 


N 


M 


N 


\  o 

1 



} 

X  X 
XXX 

o 

o 

o 

XXX   X 

xxx  • 

xxx 

XXX 

xxx 

o 

o 

xxx 

xxx 

X 

xxx 
xxx 
xxx 

xxx 
xxx 
xxx 

x  o 

X  X  w 

xxx 
xxx  x 

o 

xxx 

xxx 

xxx 

\  o 

xxx 

xxx 

xxx 

* 

xxx 

xxx 

• 

xxx 

o   \ 

In  practice  the  matrix  need  not  be  physically  sorted,  all  that  is 
needed  are  the  indices  of  the  pivot  positions. 

The  NN  "block  contains  logical  vectors  in  the  basis;  these  will 
produce  no  T]  records.   The  RR  and  CC  blocks  are  triangular-   The  MM  block 
contains  rows  and  columns  which  cannot  be  triangularized  by  sorting.  LU 
decomposition  is  applied  on  this  block. 

Rewrite  the  matrix,  in  partitioned  form,  as : 


" 

R 

0 

0 

S 

M 

0 

C 

N 

D 

I 
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where  (I  and  I  are  zero  and  identity  matrices  of  appropriate  sizes.   This 
decomposes  to: 


R 

0 

0 

I 

i — r— 

I 

0 

T 

0 

X 

M 

u 

0 

X 

I 

s 

0 

I 

0 

N 

I 

0 

0 

c 

0 

I 

0 

I 

D 

I 

The  first  and  last  matrices  are  lower  triangular,  pivoting  down  the  diagonal 
gives  T\   vectors  with  no  increase  in  the  number  of  nonzero  elements. 

If  M  =  LU,  L  and  U  being  lower  and  upper  triangular  respectively, 
then 


I 

0 

0 

0 

0 

M 

N 

I 

0 

I 

r— 
I 

0 

0 

0 

0 

L 

0 

I 

0 

I 

X 


I 

0 

0 

0 

0 

u 

N 

I 

0 

I 

The  first  matrix  is  lower  triangular  and  the  second  can  be  transformed  by  row 
and  column  interchanges  into: 


I 

0 

0 

0 

0 

I 

N 

0 

I 

0 

U 

37 


which  is  upper  triangular.  After  inversion,  the  interchanges  can  be  reversed. 
The  only  place  where  elements  are  added  to  the  inverse  is  in  the  formation  of 
L  and  U.   Since  the  pivoting  strategy  within  the  M  block  may  be  employed, 
this  may  be  used  to  minimize  the  number  of  nonzeros  created. 

The  exact  form  of  the  pivoting  strategy  depends  on  such  factors 
as  the  number  of  columns  of  the  matrix  that  can  be  stored  in  memory  at  one 
time.  An  important  point  is  to  chose  the  pivots  dynamically,  at  each  step 
choosing  that  pivot  which  minimizes  the  number  of  nonzeros  created. 

It  may  be  necessary  to  change  the  pivot  selection  method  where 
choosing  a  particular  pivot  would  result  in  very  large  elements  in  the  repre- 
sentation of  the  inverse.   Discussion  of  this  point  is  deferred  until  the 
next  chapter. 
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6.   ERROR  ANALYSIS  OF  THE  REINVERSION  ROUTINE 


In  this  chapter  the  computation  of  the  inverse  of  the  matrix 


B   = 


R 

0 

0 

0 

S 

M 

N 

c 

D 

I 

r  s  t 

columns   columns   columns 


using  the  product  form 


i-l 


.2  „1 


T9  T5        ...T     T      ,     p   =  r  +   s  +  t, 


is  considered.   It  is  shown  that 


2  „1 


T9  T9'      .  .  .  T  T  (B  +  6B)  =  I 


i.e.,  T^  .  .  •  T  is  the  exact  inverse  of  some  perturbed  matrix  B  +  5B,  and 
a  bound  for  j  | &b| |  is  established.  As  before  B  is  split  into  the  form 
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B  = 


1 
R 

0 

r— r- 

0 

0 

S 

I 

0 

I 

0 

1 

X 




I 

r 
0 

0 

■ 

0 

0 

M 

N 

I 

0 

I 

*R 


X 


X 


I 

0 

0 

| 

0 

0 

I 

0 

c 

D 

I 

M 


Consider  first  the  inversion  of 


B 


11 


21 


31 


rl 


22 


32 


r2 


33 


b  -■  -,        b  -,    n 
r+1  1    r+1  2 


ml       m2 


rr 


.  b 


r+1  r 


mr 


J 


by: 


The  T\   vector  for  the  j   elementary  matrix  T   ( j=l(l)  r)  is  given 


T\3.   =   fl  (1/b.  .)  = 


JJ'    b..(l+S..) 
33  33 


nJ  =  n  (-b. .  x  r\j)=    "bi.i  ^1+6i,j^ 

10    J    >a (^ , 


f or  i  >  j 


Uo 


where    I S -  - 1  <  e, > 


The  terms  b   (1+6  )  may  be  regarded  as  the  elements  of  a  matrix;  hence 
ij       ij 


where 


.  .  T  is  the  exact  inverse  of  the  perturbed  matrix  (B  +  SB^) 


R- 


(*Vij 


<V« 


d  =  1(1)  r 

i  =  3(1)  m 

otherwise 


Similarly,  the  product  of  the  elementary  matrices  produced  for  the  inversion 

of  B  is  the  exact  inverse  of  (Bc  +  6BC) 

where 


(*Vij 


'5ij  <Bc>u 


j  =  r  +  s  +  l(l)p 

i  =  j(l)m 
otherwise 


and 


hi^h 


This  leaves  only  the  matrix 


BM' 


I 

0 

r~~~ 

0 

0 

0 

M 

N 

I 

0 

I 

1+1 


to  be  inverted.   Since  the  elements  of  N  are  simply  copied  into  the  7]  vectors 
with  a  change  in  sign,  the  only  error  is  that  in  inverting  M.   Therefore  set, 


M=M1  = 


in 


m 


11 

1 
21 


m 


s2 


1 
"12 


m 


22 


m 


s2 


1 

mis 


m 


2  s 


m 


ss 


12        s   s+1 
and  compute  a  sequence  of  matrices  M  ,  M  ,    .  . . ,   M  ,   M   ,  by  forming  a  new 

T\   vector  from  the  k-th  column  of  NT  using  only  elements  on  or  below  the 

diagonal. 

Thus 


*£-  (1+6^) 


kk 


Jkk; 


i 


fK-4^)=-4(i+8^) 
"4  (i+4> 


where 


Is, 


ik1 


<*{■ 


(i) 


The  matrix 


M*  i 


s  then  updated  by  this  new  elementary  matrix  giving 


k+1  v+-i 

M   .   If  there  were  no  rounding  errors,   the  elements  of  M    would  be  zero 

k+1 
in  the  k-th  column  below  the  diagonal  and  el.   Would  be  unity.   If  we  set 

these  to  zero  and  one  respectively  then: 


1+2 


i.e. 


Thus  if 


m. 


k+1 

i.  . 


0 


fi  cn. 


k  \ 


^Ki^i^) 


k 
m.  . 


i  =  j  =  k 
i  >  k,  j  =  k 

i  =  k,  j  >  k 
1  >  k,  j  >  k 

otherwise      (2) 


k+1 


m 


=< 


\  \,  +  \  \j  Bk! 

+  T)*mJL  +aJd  Bid  +  2  ^i  Vi  K* 


i  =  j  =  k 
i  >  k,  j  =  k 
i  =  k,  j  >  k 


ij 


m. 


iJ 


i  >  k,  j  >  k 

otherwise 


(3) 


k 
k 
k 
k+1 


Tl 

it 


and  ET  is  the  matrix  with  elements  6.., 


then 


Mk+1      =     lV+f? 


w 


where 


U3 


k     nk    *k 
"kk^k    Skk 

k       k      L     k    ~k     k 
mik  5ik  +  ffikkT1i  *kk 

k    _k       k 


d    k     kJ 

k       k  _     k    „k      ,k 


i  =  j   =  k 
i  >  k,    j   =  k 

i  =  k,   j  >  k 
i  >  k,  j  >  k 

otherwise 


(5) 


and 


4's* 


6i*ls4> 


for  all  i,   j,  k. 


s+1 
After  s   iterations  M         =  U,    (say),    an  upper  triangular  matrix  which 

has   ones  down  the  diagonal.      Let 


F 

^+1 


-  UV1  Ek 

-  Lk  M*  +  Lk  Fk 
=  Lk  (t^  +  Pk) 


k  k 

But  the  elements  of  F  are  zero  except  for  i  >  k  and  j  >  k  (because  those  of  E 

are) . 


L1  Fk  =  Fk   f or  i  <  k 

.-.    l*"1  Lk-2  .  .  .  L2  L1  ^  -  Fk 


kk 


then 


U  =  M 


s+1 


.  s  s-1 
=  L  L 


^L1 


M1  \ L  (L$  Is"1  . 

k=l  v 


•  L  F  ) 


LS  L3"1  .  .  .  I2  L1  M1  ^  (LS  L8"1  .  .  .  L2  L1  Fk) 


LS  L3"1  .  .  •  L2  L1  (M  +  F) 


(6) 


where 


If 


■A  **• 


U 


^2     ^3 


u 


23 

1 


•l 


U2s 


ii 


u 


3s 
s-ls 

1 


u"1  =  u2  U3  .  .  .  u8"1  us 


where 


"*5 


V1   = 


■\: 


■U2j 


-U.  -  . 

J-lJ 


from  which 


U2  U3  .  .  .  U8-1  US  LS  L5"1  .  .  .  L2  L1 


L  L  (M+F)  =  I. 


A  bound  on  | |f| |  is  now  required.   Since 


F*  -  (lk)"X  ^  , 


k   k 

m.,  8., 
lk  lk 


ij 


I 


k   k 

id  ^j 


k   k    _  „k  k  ,  ,k 
m.  .  8.  .  +  2  T] .  m,  .  (5  . 


i  >  k,   j  =  k 

i  =  k,   j  >  k 

8M> 

i  >  k,   j  >  k 

otherwise 


(7) 
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An  important  point  is  that  these  error  terms  do  not  involve  the  'pivot' 
element  T|,    Summing  over  k  gives 


F.  . 


j-l 


m 


J    v  /  k   k    „  -k  k  ,„  ,k   Ak  XN 


ij   ij 


1J     KJ 

i  >  J 


< 


i-1 


z  k   k    _  ^k  k  /  ,k    k  N  N 

mij  Bij+k§i(mij  h3  +  2\\*  Kj-  "^ 


i  <  J 


(8) 


Let 


P 

M 


a.  . 


max    {   m.  . 


max    (  |T]_.  I  ) 
i,k 


4i 


Hi     k=l     l 


Z,  V  (mk.) 


i  >  k  } 


Then 


|F.  .1  <  (cr.  .  +  kq.)   p  M  £' 


|  |F|  |  <  s  (a  +  l+q)  p  M  e| 


where     q    =  max  (q  ) 

i 


=  max  ( cr .  . ) 


1,J 


ij 


This  gives  an  a  posteriori  bound  for  | |f| |  in  terms  of  the  sparse- 
ness  and  the  size  of  the  elements  of  the  inverse  and  the  intermediate  matrices 
M1,  V? ,    •  •  •,   $ '     As  might  be  expected,  sparse  matrices  give  much  smaller 
errors. 


1*7 


Combining  these  results  indicates  that  the  computed  inverse  is  the 
exact  inverse,  in  product  form,  of  the  perturbed  matrix 


B  +  5B 


1 

R+6R 

0 

0 

0 

S+8S 

M+F 

N 

C+8C 

'  D+6D 

I 

It  is  clear  from  equation  (9)  that  both  sparseness  and  size  of  the 
elements  are  crucial  factors  in  bounding  the  error.   The  previous  chapter 
dealt  with  preserving  sparseness.   Consideration  now  must  be  given  to  bounding 
the  size  of  the  elements  in  the  T)  vectors  and  the  intermediate  matrices  that 

are  formed.   In  forming  an  T)  vector,  start  with  some  column  vector 

T 
v  =  [v^,  v  ,  ...,  v  ]   (say)  and  'pivot'  on  some  element  v  ,  giving 

m 

T[   =  fl  L-v  /v  ,  .  .  .  ,  -v   /v  ,  1/v  ,  -v  _/v  ,  • . .,  -v  /v  ] 
1'  r         r-1'  r   '  r    r+1'  r        v     r 


where  the  value  of  r  may  be  arbitrarily  chosen.   The  usual  method  is  to  select 
the  pivot  element  so  that  j v  |  is  greater  than  some  'tolerance'  value  so  that 
the  elements  of  T]  will  not  be  'too  large'.   Since,  as  equation  (7)  shows, 
the  element  l/v  does  not  appear  in  any  error  term,  the  remaining  terms  may 


U8 


be  bounded  by  c  (say)  by  insisting  that  r  be  chosen  such  that: 

|vrl  >  \     |v.|  for  1-1,  ...,  r-1, 

r+1,  ...,  p 

It  may  be  noted  that  for  dense  systems,  one  normally  selects  the  largest 
element  of  v  for  this  pivot,  which  is  equivalent  to  taking  c  =  1.   Thus  the 
above  scheme  is  a  compromise  between  minimizing  the  number  of  nonzero  elements 
and  reducing  their  size. 


h9 

7.   ITERATIVE  REFINEMENT 

Suppose  that  an  approximate  inverse  to  the  basis  B  and  an  approxi- 
mate solution  vector  x  have  been  calculated.   It  is  desired  to  calculate  a 
more  accurate  solution  using  iterative  refinement.  A  detailed  analysis  of 
iterative  refinement  using  floating  point  arithmetic  is  given  by  Moler  [7] • 
A  simple  adaption  of  that  work  is  developed  here  to  cover  the  sparse  matrix 
case. 

Each  iteration  of  the  refinement  requires  the  following  three  steps 

1)  Compute  the  residuals:   r  =  b  -  Bx  with  x  =  x. 

2)  Compute  the  correction:   d  =  B  r  . 

^\   « ^  j_!        j_-      j>  21+I    m   .m 

3)  Add  the  correction  to  form:  x    =  x  +  d  . 

The  step  where  most  accuracy  is  required  is  in  calculating  the  residuals. 
It  is  assumed  that  extended  precision  is  used  in  this  step.   Computing 
d  =  B  r  is  done  using  the  approximate  inverse  to  B. 


Let 


rm   =  b  -  B  xm  =  cm  W 

dm       =     (B  +  E111)    -1  rm  (2) 

xm+l=     /i+dm+gm  (3) 


mm  m 

where   c   ,   E     and  g     are   correction  terms  due  to  round-off  error  and  are   intro- 
duced to  make  the  equations  hold  exactly.     A  bound  for   ||E  ||    is  known,   using 
the  results   of  either  Chapter  h  or  Chapter  6.      It   is  required  that  the  inverse 
be  reasonably  good  so  that 


-1  I    I    N  "I 


A\  <  (Mb"1!!) 


Let 


Fm  =  B'1  E"1. 


Then 


|Fm||  =  Mb'1  ^\\  <  Mb"1!!  ||Em||  <i. 


and  hence  I  +  F^  is  nonsingular. 
Equation  (2)  becomes 


Suppose  x  is  the  exact  solution  x  =  B  "'  b 


Then 


m+1  m  mm 

x  -x  =  x     -x  +  d     +g 
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=   (I  +  FT1  B_1   rm.  (k) 


(I+fV   [(I+^)    (xm-x)   +  B-^b-Bx10  +  c*)]   +  Sm 


-1   rjn/   m     x  -1     mn  m 


=  d^y1  [/(/-x)  +B--Lcjn]  + 


Hence 


xm+1-x||<       11^1  I        I.  m     ,,  +     Mb"1!!      licPii  +  MJni 

*  I  I  —    K 

1-1 ll"! I  l-||Fm|| 
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Suppose  an  upper  bound  for  \\^\\    is  given  by  |  |  gB  |  |  .   Then 


^'*,lsrrj{SW(,,,B,''l'--''  +  ii-"''>  +  Mfli 


Estimates  of  | | cm| |  and  ||gm||  are 


now  required, 


m 


=  fl  (xm  +  dm)  -  (xm  +  dm) 


h  =  e±\  +  e[^    >    ^r  |e.|<ei,|e.|<ei 


ui 


Is  N  <\(UAl  +  l|dffi||) 


When  the  process   converges,    ||dm||   <   | |xm|  I 


Kg"1!!   <2ex    ||xm|| 


<2ex   (||x||    +   ||xm-x||). 


m 


In  calculating  r  ,   it  is   assumed  that  each   nnrrmnnpnt   J* 


assumed  that  each   component  r.    =  fl  (b.-Z  B. .  xm) 

is   calculated  using  extended  precision  and  the  result  then  is  rounded'to'nor- 

mal  precision.      This   is  denoted  by  fl  (b      -  Z  B       xm) 

2  A        i  i.i      .1'' 


ij      J 


m         m  m 

c     =  r     -  b  +  Bx 


and 
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m 
c.    = 


n        (b.  -Zb..  xm)  -  (b.  -Lb..  xm) 

2,1   v    i  ij      j'  i  yd 

.  (fi2  (bt  -  Z By  *»))  (i  +  ej  -(b.  -Eb..  xp 

=      (b.(l^    .)       -    Z     By      *"      (1      +     7y)l        (1      +      9^       -      (b.       -     Z     By      X^) 
.      (b,       -     Z      By       X")        ,± 


where  \e   \   <  e  ,    |£±|   <  e, 

|7.J    <  (q  +  2)   e^ 

m 

q  =  max.L.   V     (B     ) 


since 


||cm||    <€^      ||b    -  Bxm    ||    +   (    1  +€]L)   ^    ||b|| 

+  (q+  2)  (l  +  ex)  e<     ||b||     ||xm|| 

<€'     ||B||       ||xm   -   x||       +  (l  +  6X)    ^      llBH     l'Xl 

+  (q+  2)  (l  +  ex)  e2    ||b||  (||x||  +  ||xm  -  x||) 

i  i    nil  i         M  /   m  \  i  i 

| |x    | |    =    | |x  +   (x     -  x)| I 

I  I    I  I        I  I   m         II 
<        |x| I    +    I |x      -  x| I . 

I  <    [€'  +  (q  +  2)  (i  +  ex)  €']    ||B|| 
+  (q+  3)   (1  +  ex)  €2     HB"    "X"' 


m  i 

x     -  x 


m+1 


-  x|  I    <    I  |5B| 


-1 


m  i 

x     -  x 


1    -    |  l&B 

+  1 1^11     INI     M^M 
l  -  I  IfiP'l  I     I  |b"  1 1 


(e£   +    (q  +    2)    (1   +  €1)    G2) 

+  (q+  3)   (i  +  e{)  €'    I 


m  I 

x      -  x 


+   2 


ex  (IM 


m  I 

x     -  x 
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SB | |     | |B      | |    +  k  (B)    (6'    +   (q  +  2)    ( 


1  -       SB 


-1] 


1  +  g1)  e£)    +  2<H 


xm-xl 


JlM 


1    -    ,    5B  B 


—    (q+  3)   (1  +  e1)  e£  +  2e1 


] 


where     k  (B)   =   |  |b|  j     ||b  ~  |  |  ,    the   condition  number  of  B. 


I  I    m  II  Mil  /  \ 

0l   |  |x     -  x|  I    +   a2\  lxl  I        (say). 


m+1        i|.  if   m       ii 

x         -x  |  |     <     cr       j  I  x     -x  I  j      + 


a2 


x 


Since  x     is   calculated  in  the   same  way  (using,    in  effect, 


x°   =     0) 


x     -  x| I   <     CTl    I |x| 


(This  may  be  shown  by  setting  e'  =  0  for  this  step  (no  calculation  of  re- 
siduals is  involved)  and  noting  that  the  only  other  term  in  a2t   2e  f    came  from 
2ei    ||xm||  <  2€  (||x||  +  ||xm  -  x||),  but  x°  =  0). 


1    I  i 
x  -  x|  1  <  cr 


x 


by  induction, 


I    m  I 

x     -  x 


m 


<      ul 


On        +  0^  >  O,       <     1 


1  - 


1_0H 


5h 


Therefore,  provided  a <   1,  convergence  of  x  to  within  the 


tolerance  given  by 


*     11-11. 


l-°i 


can  be  guaranteed. 

cr  depends  essentially  on  |  |$b||,  lB  and  the  measure  °f  sparse- 
ness,  q.  It  is  essential  that  B  be  well  conditioned  and  that  ||6B||  be  small, 
i  e.,    an  accurate  inverse  of  B  is  known. 

The  effect  of  a  small  value  of  q  is  to  reduce  the  need  for  using 
double  precision  in  calculating  the  residuals,  since  it  is  always  multiplied 

by  ^ 


8.   CONCLUDING  REMARKS 

The  question  remains  as  to  what  constitutes  an  optimal  feasible 
solution.  Strictly  speaking,  it  is  one  for  which  x  >  0  and  the  reduced  costs 
d  are  also  >  o  for  j  =  l(l)n.   In  practice,  small  negative  tolerances,  5  ,  5 
(say)  are  chosen  so  that  the  solution  is  considered  optimal  and  feasible  pro- 
vided 


x.  >  g  for  all  i 


dj  >  SQ  for  all  columns  j  not  in  the  basis 


An  alternative  approach  is  as  follows:   Determine  some  sufficiently 
refined  solution  x  and  a  basis  B,  such  that 


Bx  =  b 


may  be  considered  exact.   If  some  x.  are  <  0,  define 

1        ' 


x  =  x  +  6x  >  0 


where 


Bx. 


CO        ,   x.> 
\Xi  '  Xi< 


0 


This  produces  a  change  in  b,  5b  (say),  where 
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6b  =  B  6x 

The  vector  5b  can  either  be  calculated  directly,  or  its  norm  may  be  estimated 
from 

||sb||  <  ||B||  ||fix|| 

If  the  elements  of  8b  are  sufficiently  small  that  b  +  5b  is  an  acceptable 
right  hand  side,  then  x  may  be  considered  to  be  the  required  solution  to  the 
problem. 

Similarly,  in  calculating  d.  =  c.  -.A.  n.  a.  .,  n  may  be  refined  in 

the  same  way  as  x,  so  that  the  d.  are  accurately  known.   If  d.  <  0  for  some  j, 

J  J 

a  small  change  in  the  corresponding  c.  will  set  them  to  zero,   if  this  new 

J 

cost  vector  is  acceptable,  then  the  problem  is  solved. 

The  techniques  described  in  Chapters  3  through  6  and  above  enable 
the  exact  solution  to  the  slightly  perturbed  problem; 

minimize 

T 
z  +  5z  =  ( c  +  6c)  x 

subject  to 

(A  +  gA)  x  =  b  +  6b 
x  >  0. 

to  be  obtained.   If   |  oA|  |  ,   |  &t> |   and  ||5c||  are  within  acceptable  limits, 
then  the  problem  may  be  considered  solved.   If  more  accuracy  is  required, 
iterative  refinement  as  described  in  Chapter  7  may  be  used. 


57 
LIST  OF  REFERENCES 


[1]  Bartels,  Richard  H. ,  "A  Numerical  Investigation  of  the  Simplex  Method", 
Technical  report  no.  CS  104,  Stanford  University,  1968. 

[2]   Clasen,  R.  J. ,  "Techniques  for  Automatic  Tolerance  Control  in  Linear 
Programming",  Comm.  ACM,  £  (November  I966),  pp.  802-803. 

[3]  Dantzig,  George  B. ,  Linear  Programming  and  Extensions,  Princeton 
University  Press,  I965. 

[U]  Dantzig,  George  B.,  Harvey,  Roy  P.,  McKnight,  Robert  D. ,  Smith,  Stanley  S. , 
"Sparse  Matrix  Techniques  in  two  Mathematical  Programming  Codes", 
Technical  report  no.  69-I,  Stanford  University,  1969. 

[5]  Forsythe,  George  E.  and  Moler,  Cleve  B. ,  Computer  Solution  of  Linear 
Algebraic  Systems,  Prentice  Hall,  I967. 

[6]  Hadley,  G. ,  Linear  Programming,  Addison  Wesley,  I962. 

[?]  Moler,  Cleve  B. ,  "iterative  Refinement  in  Floating  Point",  JACM,  1^ 
(April  I967),  pp.  316-321. 

[8]   Orchard-Hays,  William,  Advanced  Linear  Programming  Techniques,  McGraw 
Hill,  1968. 

[9]  Wilkinson,  J.  H. ,  Rounding  Errors  in  Algebraic  Processes,  Prentice  Hall, 
1963. 


58 


APPENDIX 
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Lemma  1;  The  solution  of  a  triangular  system  of  equations. 

Suppose  that  the  components  of  a  vector  x  are  computed  in  the  order 
\>  x2f    '"    '   Xn  by 


x,  =  b 


1    1 


x.  = 


i-1 

fl  (b.  +  JL    rM  x,  ) 
i   k=l  ik  k' 


i =  2,  3,    •••  ,  n 


or,  equivalently,  by  the  set  of  equations 


■r21    1 


■r     -r 
31    32 


-r     -r 
ml     m2 


-r 


ram-1  1 


x~ 


LXm- 


m 


i.e.,  Rx  =  b,  is  solved.   Then  x  is  the  exact  solution  of  the  perturbed  tri- 
angular system: 

(R  +  SR)  x  =  b  +  sb, 


where, 


6R  = 


-r21821        0 
'r31531    "r32632 


ml  ml 


0. 


mm-1  mm-1 
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and 


where 


6b  =  [bx  ai ,  •  •  • ,  \  s;]T , 


Proof 


s:  I  <  «#      - 


i-l 
with  ^  ■&  V     (rik) 


In  Chapter  2,  it  was  shown  that 


n 


fl  &  vi}  -A  v±  <1+8i> 


Bj.  I  <  (^c  +  1)  €' 


n 


and  c^  -J^  V  (x±) 


Similarly,  if  p  =  f  1  ( z  +2^  x±y±)    is  calculated  from: 


Sl  =  Z 


B.  .  =  fl  (s.  +  xiyi)  i  =  1,  2,  ...  ,  n 


p  =  Sn+1 


then 


n  n 

fl  (z  +i§i  xiyi}  =  2  (1  +  ^)  +jjl  x.y.  (1  4  6.) 


where 


A  |  <  q^  e«   , 


I  $!    I  <  (%,  +   1)  €' 


and  qx  =.Z1  V  (x.)   . 


Using  this  result, 


Xi   ■  fl   (\  +^l      riK  *J 


"  ^   (1  +   ^      +  li   **  **  (1  +   ^ 


where 


Bik   I   <(<11  +  1)€^ 


and 


8.    |   <  q.    el      • 


6l 
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Hence 


i-1 
Xi  =  bi  +  \  6i  \ll   (rik  +  rik  &ik>  \ 


which  completes  the  proof. 
Note  also  that 


i-1 
||5R||  ^max.^  |  r.  .  8.. 


<  max  (q         max        |r     8     |  ) 
i         x  Kj<i-1       1J   1J 


<  q  (q+l)   €'   max      |r. . 


where   q  =  max  q. 
i        1 


If  r   results  from  a  division,  as  described  earlier,  two  extra  errors  are 

ik 

involved  so  that,  in  the  lemma, 


6ikl  <  (^  +  3)  «i 


and  |oR||    <     q  (q  +   3)   e'      max        |r 
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Lemma  2.  A  direct  error  bound  on  x. 

Under  the  same  hypotheses  as  in  Lemma  1,    if  x  is  the  exact  solution 
R  b,  and  8*  =  x  -  x,  then 


1 6x1  I  <    l|R  "    [m||x||  +  e«  ||b||] 


where         m  =  q  (q  +  2)  e'  max   Jr.,  I   , 


i,k    ik 


provided  that 


m 


llR"-1!!  <  i    • 


Proof 


As  before, 

(R+5R)   x  =  b  +   5b 

also  R  x  =  b 

b  +  6b  =  (R  +   6R)    (x  -   5x) 

=  Rx  +  6Rx    -   (R  +  8R)    5x 
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6b   =   8Rx   -   (R  +   gR)    6x 


8x  =  (R  +  5R)"1   [6Rx  -  0b] 


-1 


and     ||&x||   <   ||(   R  +  5R)"    ||    [||SR||    |  |x|  |    +    |  |  8b|  |  ] 


but  ||(R  +  8R)"1!!    =   ||[R(I  +  R^SR)]"1 


<   ||R_1||    ||(I  +R_1   gR)-1!! 


R-\ 


l-llR^SRl 


IR"1! 


< 


i-llR"1!!   IISRl 


provided  | |r-1| |    | | &R | |    <  1 

( and  hence  | |r     &r| |<  l) 

•  ||5x||    <     [J 


R"1! 


,— —  UMI    ||x||    +  e«    ||b||] 

i-IIr    II  IIsrII  1 


But  x  =  x  -  5Xj 


x   :   x   +   5x 


|ox||.   N^ 


,,  -lM  ,,   ,,  Ell6R||  (||x||  +  ||*||)  +€J  ||b||] 
1-||R  II   6R  X 
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which  reduces  to 


6x11  <  U**1 


,,  -i,,  ,,    ,,    rllsBll  11*11  +  *i  IMIi 
i-2||k    II  IIsrII  l 


provided 


■"ll    IISRll   <i 


and  m  =   |  |  SR  | 
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